function qr_experiment()
% QR_EXPERIMENT QR分解综合实验
% 
% 本实验涵盖QR分解的各种算法和最小二乘应用

fprintf('=== QR分解综合实验 ===\n\n');

%% 实验1：Gram-Schmidt QR分解
fprintf('实验1：Gram-Schmidt QR分解\n');
fprintf('目标：理解QR分解的基本原理\n\n');

% 测试矩阵
A1 = [1, 1, 0; 1, 0, 1; 0, 1, 1; 0, 0, 1];
fprintf('测试矩阵 A1 (%d×%d) =\n', size(A1));
disp(A1);

% Classical Gram-Schmidt
[Q1_cgs, R1_cgs] = gram_schmidt_classical(A1);
fprintf('Classical Gram-Schmidt结果:\n');
fprintf('Q矩阵 =\n'); disp(Q1_cgs);
fprintf('R矩阵 =\n'); disp(R1_cgs);

% 验证正交性
orthogonality_cgs = norm(Q1_cgs' * Q1_cgs - eye(size(Q1_cgs, 2)), 'fro');
fprintf('Q的正交性误差: %e\n', orthogonality_cgs);

% 验证分解
reconstruction_cgs = norm(Q1_cgs * R1_cgs - A1, 'fro');
fprintf('重构误差: %e\n', reconstruction_cgs);

% Modified Gram-Schmidt
[Q1_mgs, R1_mgs] = gram_schmidt_modified(A1);
fprintf('\nModified Gram-Schmidt结果:\n');
orthogonality_mgs = norm(Q1_mgs' * Q1_mgs - eye(size(Q1_mgs, 2)), 'fro');
reconstruction_mgs = norm(Q1_mgs * R1_mgs - A1, 'fro');
fprintf('Q的正交性误差: %e\n', orthogonality_mgs);
fprintf('重构误差: %e\n', reconstruction_mgs);

%% 实验2：数值稳定性比较
fprintf('\n实验2：数值稳定性比较\n');
fprintf('目标：比较Classical和Modified Gram-Schmidt的稳定性\n\n');

% 构造病态矩阵
epsilon = 1e-7;
A2 = [1, 1; epsilon, 0; 0, epsilon];
fprintf('病态矩阵 A2 (ε = %e):\n', epsilon);
disp(A2);
fprintf('条件数: %.2e\n', cond(A2));

% Classical Gram-Schmidt
[Q2_cgs, R2_cgs] = gram_schmidt_classical(A2);
orthogonality2_cgs = norm(Q2_cgs' * Q2_cgs - eye(size(Q2_cgs, 2)), 'fro');

% Modified Gram-Schmidt
[Q2_mgs, R2_mgs] = gram_schmidt_modified(A2);
orthogonality2_mgs = norm(Q2_mgs' * Q2_mgs - eye(size(Q2_mgs, 2)), 'fro');

fprintf('病态矩阵的正交性误差:\n');
fprintf('Classical GS: %e\n', orthogonality2_cgs);
fprintf('Modified GS: %e\n', orthogonality2_mgs);
fprintf('改善倍数: %.1f\n', orthogonality2_cgs / orthogonality2_mgs);

%% 实验3：Householder QR分解
fprintf('\n实验3：Householder QR分解\n');
fprintf('目标：掌握最稳定的QR分解方法\n\n');

A3 = randn(6, 4);
fprintf('随机矩阵 A3 (%d×%d)\n', size(A3));

% Householder QR分解
[Q3_house, R3_house] = householder_qr_simple(A3);
fprintf('Householder QR分解完成\n');

% 验证结果
orthogonality3 = norm(Q3_house' * Q3_house - eye(size(Q3_house, 2)), 'fro');
reconstruction3 = norm(Q3_house * R3_house - A3, 'fro');
fprintf('Q的正交性误差: %e\n', orthogonality3);
fprintf('重构误差: %e\n', reconstruction3);

% 与MATLAB内置QR比较
[Q3_matlab, R3_matlab] = qr(A3, 0);
Q_diff = norm(abs(Q3_house) - abs(Q3_matlab), 'fro');  % 符号可能不同
R_diff = norm(abs(R3_house) - abs(R3_matlab), 'fro');
fprintf('与MATLAB QR的差异: Q差异 = %e, R差异 = %e\n', Q_diff, R_diff);

%% 实验4：最小二乘问题求解
fprintf('\n实验4：最小二乘问题求解\n');
fprintf('目标：使用QR分解求解最小二乘问题\n\n');

% 多项式拟合问题
x_data = linspace(0, 1, 10)';
y_true = 2*x_data.^2 + 3*x_data + 1;
noise = 0.1 * randn(size(x_data));
y_data = y_true + noise;

% 构造Vandermonde矩阵（2次多项式）
A4 = [ones(length(x_data), 1), x_data, x_data.^2];
fprintf('多项式拟合问题:\n');
fprintf('数据点数: %d\n', length(x_data));
fprintf('多项式次数: 2\n');
fprintf('系数矩阵条件数: %.2e\n', cond(A4));

% 使用QR分解求解
[Q4, R4] = qr(A4, 0);
coeffs_qr = R4 \ (Q4' * y_data);
fprintf('QR分解解: [%.4f, %.4f, %.4f]\n', coeffs_qr);
fprintf('真实系数: [1.0000, 3.0000, 2.0000]\n');

% 与正规方程比较
coeffs_normal = (A4' * A4) \ (A4' * y_data);
fprintf('正规方程解: [%.4f, %.4f, %.4f]\n', coeffs_normal);

% 残差分析
residual_qr = norm(A4 * coeffs_qr - y_data);
residual_normal = norm(A4 * coeffs_normal - y_data);
fprintf('残差比较: QR = %.6f, 正规方程 = %.6f\n', residual_qr, residual_normal);

%% 实验5：性能基准测试
fprintf('\n实验5：性能基准测试\n');
fprintf('目标：比较不同QR算法的性能\n\n');

sizes = [100, 200, 500];
fprintf('矩阵规模    MGS        Householder    MATLAB\n');
fprintf('--------    ---        -----------    ------\n');

for i = 1:length(sizes)
    m = sizes(i);
    n = min(m, 50);
    A_bench = randn(m, n);
    
    % Modified Gram-Schmidt
    tic;
    [Q_mgs, R_mgs] = gram_schmidt_modified(A_bench);
    time_mgs = toc;
    
    % Householder QR
    tic;
    [Q_house, R_house] = householder_qr_simple(A_bench);
    time_house = toc;
    
    % MATLAB内置QR
    tic;
    [Q_matlab, R_matlab] = qr(A_bench, 0);
    time_matlab = toc;
    
    fprintf('%6d      %6.4f     %10.4f     %6.4f\n', ...
            m, time_mgs, time_house, time_matlab);
end

%% 实验6：应用实例 - 数据拟合
fprintf('\n实验6：应用实例 - 数据拟合\n');
fprintf('目标：使用QR分解进行曲线拟合\n\n');

% 生成实验数据
t = linspace(0, 2*pi, 20)';
y_exact = sin(t) + 0.5*cos(2*t);
y_noisy = y_exact + 0.2*randn(size(t));

% 构造傅里叶基函数
n_terms = 5;
A_fourier = ones(length(t), 1);
for k = 1:n_terms
    A_fourier = [A_fourier, cos(k*t), sin(k*t)];
end

fprintf('傅里叶级数拟合:\n');
fprintf('数据点数: %d\n', length(t));
fprintf('基函数数: %d\n', size(A_fourier, 2));

% QR分解求解
[Q_fourier, R_fourier] = qr(A_fourier, 0);
coeffs_fourier = R_fourier \ (Q_fourier' * y_noisy);

% 计算拟合结果
y_fitted = A_fourier * coeffs_fourier;
mse = mean((y_fitted - y_exact).^2);
fprintf('均方误差: %.6f\n', mse);

% 绘制结果
figure('Name', 'QR分解数据拟合');
plot(t, y_exact, 'g-', 'LineWidth', 2);
hold on;
plot(t, y_noisy, 'ro', 'MarkerSize', 6);
plot(t, y_fitted, 'b--', 'LineWidth', 2);
xlabel('t');
ylabel('y');
title('基于QR分解的傅里叶级数拟合');
legend('真实函数', '噪声数据', 'QR拟合', 'Location', 'best');
grid on;

%% 实验总结
fprintf('\n=== 实验总结 ===\n');
fprintf('1. Modified Gram-Schmidt比Classical更稳定\n');
fprintf('2. Householder QR是最稳定的方法\n');
fprintf('3. QR分解在最小二乘问题中表现优异\n');
fprintf('4. 不同算法适用于不同的应用场景\n');
fprintf('5. 数值稳定性在病态问题中尤为重要\n\n');

end

% 辅助函数
function [Q, R] = gram_schmidt_classical(A)
% Classical Gram-Schmidt正交化
[m, n] = size(A);
Q = zeros(m, n);
R = zeros(n, n);

for j = 1:n
    v = A(:, j);
    for i = 1:j-1
        R(i, j) = Q(:, i)' * A(:, j);
        v = v - R(i, j) * Q(:, i);
    end
    R(j, j) = norm(v);
    Q(:, j) = v / R(j, j);
end
end

function [Q, R] = gram_schmidt_modified(A)
% Modified Gram-Schmidt正交化
[m, n] = size(A);
Q = zeros(m, n);
R = zeros(n, n);
V = A;

for j = 1:n
    R(j, j) = norm(V(:, j));
    Q(:, j) = V(:, j) / R(j, j);
    for k = j+1:n
        R(j, k) = Q(:, j)' * V(:, k);
        V(:, k) = V(:, k) - R(j, k) * Q(:, j);
    end
end
end

function [Q, R] = householder_qr_simple(A)
% 简化的Householder QR分解
[m, n] = size(A);
Q = eye(m);
R = A;

for k = 1:min(m-1, n)
    x = R(k:m, k);
    
    % 计算Householder向量
    v = x;
    v(1) = v(1) + sign(x(1)) * norm(x);
    v = v / norm(v);
    
    % 应用Householder变换
    H = eye(m-k+1) - 2 * (v * v');
    R(k:m, k:n) = H * R(k:m, k:n);
    Q(k:m, :) = H * Q(k:m, :);
end

Q = Q';
end